subroutine CalcViscousFluxesForCell(iCell)
use Data_Mesh, only: meshCellPoints, meshCells, meshFaceOwner, meshFaces, meshSf
use DataPerfectGas, only: nu
use DataFaces, only: FaceVarsNew, FaceVarsOld, ViscousFluxes
use Data_Cfoam, only: ZeroEpsilon
implicit none
    integer,intent(in) :: iCell

    real(8) :: GradU(3,3,8), Stress(3,3,6), ViscousFluxesLocal(3), Rho, sign
    integer :: iPointCell, iPointFace, iPoint, i,j, iFaceCell, iFace, iSide
    
    if(nu.lt.ZeroEpsilon) return
    
    ! Calculation of GradU for every nodes of cell iCell
    do iPointCell = 1,8
        iPoint = meshCellPoints(iPointCell,iCell)
        call CalcGradUforPoint(iPoint,iCell, GradU(1:3,1:3,iPointCell))
    end do

    do iPointCell = 1,8
		do i=1,3
			do j=1,3
				gradU(i,j,iPointCell) = nu*(gradU(i,j,iPointCell) + gradU(j,i,iPointCell)) ! NOTE: include only the first viscosity, VSemiletov
					! [gradU(i,j,k)] = m2/s / s = m2/s2
			end do
		end do
    end do

    ! Calculation of viscous fluxes for every faces of cell iCell
    Stress(1:3,1:3,1:6) = 0d0
    do iFaceCell = 1,6
        iFace = meshCells(iFaceCell,iCell)
        if(iCell.eq.meshFaceOwner(iFace)) then
            iSide = 1
            sign = 1d0
        else
            iSide = 2
            sign = -1d0
        end if
        do iPointFace = 1,4
            iPoint = meshFaces(iPointFace,iFace)
            iPointCell = iPoint2PointCell(iPoint,iCell)
            Rho = 0.5d0 * (FaceVarsNew(1,iSide,iFace) + FaceVarsOld(1,iSide,iFace))
            Stress(1:3,1:3,iFaceCell) = Stress(1:3,1:3,iFaceCell) + &
                                        Rho * gradU(1:3,1:3,iPointCell) / 4.0d0
        end do

        ViscousFluxesLocal = 0d0
        do i = 1,3
            ViscousFluxesLocal(1:3) = ViscousFluxesLocal(1:3) + meshSf(i,iFace) * Stress(1:3,i,iFaceCell)
        end do
        
        ViscousFluxes(1:3,iSide,iFace) = sign * ViscousFluxesLocal(1:3)
    end do

contains

real(8) function iPoint2PointCell(iPoint,iCell)
implicit none
    integer,intent(in) :: iPoint,iCell

    integer :: iPointJ

    do iPointCell = 1,8
        iPointJ = meshCellPoints(iPointCell,iCell)
        if(iPoint.eq.iPointJ) then
            iPoint2PointCell = iPointCell
            return !iPoint2PointCell
        end if
    end do

    Stop '@@Fatal Error: iPoint2PointCell has not been finished successfully!'

end function iPoint2PointCell

end subroutine CalcViscousFluxesForCell

subroutine CalcGradUforPoint(iPoint,iCell, GradU)
use Data_Mesh, only: FacesInPointsMax, Points2Faces, meshFaceOwner, meshFaceNeighbour, meshSf, &
                     meshVPoints, nInternalFaces
use DataCells, only: SubStepCellVars
use DataFaces, only: FaceVarsNew, FaceVarsOld
use Data_Cfoam, only: Unew
use dataTime, only: TimeCell, dt2Cell
implicit none
    integer,intent(in) :: iPoint, iCell
    real(8),intent(out) :: GradU(3,3)
    
    real(8) :: LinearInterpolation
    
    integer :: iFaceNum, iFace, i,j, iCellOwner, iCellNeighbour
    real(8) :: Uowner, Uneighbour, deltaU
    
    GradU(1:3,1:3) = 0d0
    
    do iFaceNum = 1,FacesInPointsMax
        iFace = Points2Faces(iFaceNum,iPoint)
        if(iFace.eq.0) exit
        
        do j = 1,3
            iCellOwner = meshFaceOwner(iFace)
            Uowner = LinearInterpolation(SubStepCellVars(1+j,iCellOwner), &
                                         Unew(j,iCellOwner),TimeCell(iCellOwner),dt2Cell(iCellOwner), &
                                         TimeCell(iCell)+dt2Cell(iCell))
            if(iFace.le.nInternalFaces) then
                iCellNeighbour = meshFaceNeighbour(iFace)
                Uneighbour = LinearInterpolation(SubStepCellVars(1+j,iCellNeighbour), &
                                         Unew(j,iCellNeighbour),TimeCell(iCellNeighbour), &
                                         dt2Cell(iCellNeighbour), &
                                         TimeCell(iCell)+dt2Cell(iCell))
            else
                Uneighbour = 0.5d0 * (FaceVarsNew(1+j,1,iFace) + FaceVarsOld(1+j,1,iFace))
            end if
            deltaU = Uowner - Uneighbour
            GradU(1:3,j) = GradU(1:3,j) + meshSf(1:3,iFace)*deltaU
        end do
    end do
    
    GradU(1:3,1:3) = GradU(1:3,1:3) / meshVPoints(iPoint)

end subroutine CalcGradUforPoint

! Linear Interpolation for (t1,t1+dt) into t
real(8) function LinearInterpolation(U1,U2,t1,dt,t)
use Data_Cfoam, only: ZeroEpsilon
implicit none
    real(8),intent(in) :: U1,U2,t1,dt,t
    
    real(8) :: alfa
    
    if(dt.lt.ZeroEpsilon)then
        Stop '@@Fatal Error: dt = 0'
    end if
    
    alfa = min(dt,max(0d0,t1-t)) / dt

    LinearInterpolation = alfa * U1 + (1d0 - alfa) * U2

end function LinearInterpolation

!subroutine CalcGradU()
!use DataCells, only: SubStepCellVars
!use Data_cfoam
!use Data_Mesh, only: nFaces, nInternalFaces, meshFaceOwner,meshFaceNeighbour, &
!                     opposingOwnerFaceLabel,opposingNeighbourFaceLabel, &
!                     meshCells,meshFaces,meshPoints
!implicit none
!
!    integer :: iFace, iBoundField
!    integer :: i,j,k, n, a,b, iPoint1,iPoint2
!
!!    ttime=ttime+dt2
!!	pNew=0.0d0
!!	Tnew=0.0d0
!!	Unew=0.0d0
!!	Enew=0.0d0
!!	rhonew=0.0d0
!! viscosity 
!
!!=====================================================
!!	gradU=0.0d0
!!	do iFace = 1,nInternalFaces ! k=
!!		do i=1,3
!!!			delta=U(i,meshFaceOwner(k))-U(i,meshFaceNeighbour(k))
!!			delta = SubStepCellVars(i+1,meshFaceOwner(iFace)) - SubStepCellVars(i+1,meshFaceNeighbour(iFace))
!!			do j=1,3
!!				vtmp = meshSf(j,iFace)*delta
!!				do l=1,4
!!					gradU(i,j,meshFaces(l,iFace))=gradU(i,j,meshFaces(l,iFace))+vtmp ! m2 * m/s
!!				end do
!!			end do
!!		end do
!!	end do
!!
!!	do iBoundaryField = 1,nBoundaryFields ! j=
!!!		if (btype(j).ne.3.and.btype(j).ne.4.and.btype(j).ne.5) then ! INLET, OUTLET, PROC
!!		if (btype(iBoundaryField).ne.5) then ! INLET, OUTLET, PROC ! btype(j).ne.4.and.
!!!			if(btype(j).eq.1.or.btype(j).eq.2)then ! INLET, OUTLET
!!			if(btype(iBoundaryField).eq.1 &
!!			        .or.btype(iBoundaryField).eq.2 &
!!			        .or.btype(iBoundaryField).eq.3 &
!!			        .or.btype(iBoundaryField).eq.4)then ! INLET, OUTLET
!!				factor=0.0d0 ! NOTE: factor=1.0d0 give mistake: wrong velocity on the boundary, VSemiletov
!!			else
!!				factor=1.0d0	! PROC boundary
!!			end if
!!			do iFace = boundaryMeshStart(iBoundaryField)+1,boundaryMeshStart(iBoundaryField)+boundaryMeshSize(iBoundaryField)
!!				do i=1,3
!!!					delta=factor*U(i,meshFaceOwner(k)) + (1.0d0-factor)*(U(i,meshFaceOwner(k)) - us(i,k)) ! Corrected BC, SK
!!					delta = factor * SubStepCellVars(i+1,meshFaceOwner(iFace)) + &
!!					    (1.0d0-factor)*(SubStepCellVars(i+1,meshFaceOwner(iFace)) - FaceVarsNew(i+1,1,iFace)) ! Corrected BC, SK
!!!					delta=factor*U(i,meshFaceOwner(k))
!!					do m=1,3
!!						vtmp = meshSf(m,iFace) * delta
!!						do l=1,4
!!							gradU(i,m,meshFaces(l,iFace)) = gradU(i,m,meshFaces(l,iFace)) + vtmp ! m2 * m/s
!!						end do
!!					end do
!!				end do
!!			end do
!!		end if
!!
!!!		if (btype(j).eq.4) then
!!!			do k=boundaryMeshStart(j)+1,boundaryMeshStart(j)+boundaryMeshSize(j)
!!!				i1=meshFaceOwner(k)	
!!!				!	vtmp=dsqrt(U(3,i1)**2+U(3,i1)**2+U(3,i1)**2)
!!!				!	localRe=Db(k)/2d0*vtmp/nu
!!!
!!!
!!!				!	factor=rwall
!!!				!*defFactor(localRe,vtmp,nu)
!!!				do i=1,3
!!!					! SK MODIFICATION FOR THE SECOND-ORDER VISCOUS WALL TREATMENT  
!!!
!!!						delta=factor*U(i,meshFaceOwner(k)) + 2.0 * (1.0d0-factor)*(U(i,meshFaceOwner(k)) - us(i,k))  ! Corrected BC, SK
!!!!						delta=6.0d0*U(i,i1)-2.0d0*us(i,opposingOwnerFaceLabel(k))
!!!
!!!					!  ANOTHER POSSIBLE WAY(1):   =4.0d0*U(i,i1)-1.0d0*us(i,opposingOwnerFaceLabel(i))
!!!
!!!					! ORIGINAL MA'S FIRST ORDER: =factor*U(i,i1)
!!!
!!!					! END SK MODIFICATION 
!!!					do m=1,3
!!!						vtmp=meshSf(m,k)*delta
!!!						do l=1,4
!!!							gradU(i,m,meshFaces(l,k))=gradU(i,m,meshFaces(l,k))+vtmp ! m2 * m/s
!!!						end do
!!!					end do
!!!				end do
!!!			end do
!!!		end if
!!
!!!	end do
!!!!   periodic
!!!	do j=1,nBoundaryFields
!!        if (btype(iBoundaryField).eq.5) then ! j=
!!			ndeltap=boundaryMeshSize(iBoundaryField)/2
!! 			do k=1+boundaryMeshStart(iBoundaryField),boundaryMeshStart(iBoundaryField)+ndeltap
!!				do i=1,3
!!!					delta=U(i,meshFaceOwner(k))-U(i,meshFaceOwner(k+ndeltap))
!!					delta = SubStepCellVars(i+1,meshFaceOwner(iFace)) - SubStepCellVars(i+1,meshFaceNeighbour(iFace))
!!					do j1=1,3
!!						vtmp=meshSf(j1,iFace)*delta
!!						do l=1,4
!!							gradU(i,j1,meshFaces(l,iFace))=gradU(i,j1,meshFaces(l,iFace))+vtmp ! m2 * m/s
!!							gradU(i,j1,meshFaces(l,iFace+ndeltap))=gradU(i,j1,meshFaces(l,iFace+ndeltap))-vtmp
!!						end do
!!					end do
!!				end do
!!			end do
!!		end if
!!	end do
!!
!!	include 'netPointsInterface.h'
!!
!!	do iPoint = 1,nPoints ! k =
!!		do i=1,3
!!			do j=1,3
!!				grad(i,j)=gradU(i,j,iPoint) / meshVPoints(iPoint) !  ! 1 / s
!!			end do
!!		end do
!!		do i=1,3
!!			do j=1,3
!!				gradU(i,j,iPoint) = nu*(grad(i,j)+grad(j,i)) ! NOTE: include only the first viscosity, VSemiletov
!!					! [gradU(i,j,k)] = m2/s / s = m2/s2
!!			end do
!!		end do
!!	end do
!!! viscosity 
!!
!end subroutine CalcGradU
